% fig6b_improved.m
clear; close all; clc;

%% 1. 模型参数
alpha = 5;
sigma = 0.8;
mu    = 0.1;
epsilon = 0.1;

% 为充分脱离瞬态，建议使用较长的迭代次数
N_total = 1000;  % 总迭代次数
N_cut   = 100;   % 去除瞬态

% 参数扫描范围：phi0 在 [-4, 6]，k 在 [0, 0.4]
phi_min = -4; phi_max = 6; Nphi = 101;
k_min   = 0;  k_max   = 0.4; Nk = 101;

phi_values = linspace(phi_min, phi_max, Nphi);
k_values   = linspace(k_min, k_max, Nk);

%% 2. 预分配矩阵存储最大 Lyapunov 指数
LE1_map = zeros(Nphi, Nk);
% LE2_map = zeros(Nphi, Nk);

%% 3. 双参数扫描，计算 Lyapunov 指数
% 注意：这里扫描时对 phi0 的循环为外层，k 的循环为内层
parfor i_phi = 1:Nphi
    phi0_now = phi_values(i_phi);
    tmp_LE1 = zeros(1, Nk);
    % tmp_LE2 = zeros(1, Nk);
    for i_k = 1:Nk
        k_now = k_values(i_k);
        % 调用 Rulkov 模型（初始状态在 Rulkov.m 中固定为 (0.01, 0, -0.01, 0)）
        [X1_save, ~, ~] = Rulkov(alpha, sigma, mu, k_now, phi0_now, N_total, epsilon, N_cut);
        % 取去除瞬态后的有效数据
        idx = (N_cut+1):N_total;
        X1_eff = X1_save(idx);
        % 调用 LEs.m 计算最大和次大 Lyapunov 指数（此处只用 LE1）
        [LE1, ~] =LEs(alpha, sigma, mu, k_now, phi0_now, epsilon);
        tmp_LE1(i_k) = LE1;
        % tmp_LE2(i_k) = LE2;
    end
    LE1_map(i_phi, :) = tmp_LE1;
    % LE2_map(i_phi, :) = tmp_LE2;
end

%% 4. 绘制最大 Lyapunov 指数谱图
figure;
imagesc(k_values, phi_values, LE1_map);
axis xy;  % 使得 y 轴从下向上递增
xlabel('Coupling strength k','FontName','Times New Roman','FontSize',12);
ylabel('Memristor initial state \phi_0','FontName','Times New Roman','FontSize',12);
title('(b) Maximum Lyapunov Exponent (LLE)','FontName','SimHei','FontSize',12);
colorbar;
colormap jet;
ylim([phi_min, phi_max]);
